##########################################################################################
##########################################################################################

library(dplyr)
library(data.table)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--sample_file"), type = "character") ,
    make_option(c("--gistic_file"), type = "character") ,
    make_option(c("--burden_file"), type = "character") ,
    make_option(c("--out_file"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    sample_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")
    burden_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.tsv",sep="")
    gistic_file <- paste(work_dir,"/gistic/all_lesions.conf_99.txt",sep="")
    ######
    out_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_file <- opt$sample_file
gistic_file <- opt$gistic_file
burden_file <- opt$burden_file
out_file <- opt$out_file

###########################################################################################

info_mss <- data.frame(fread(sample_file))
seg <- data.frame(fread(gistic_file))
dat_burden <- data.frame(fread(burden_file))

###########################################################################################

info <- subset( info_mss , Class != "IM" )
info$t_n <- paste0( info$Tumor , "_" , info$Normal )

###########################################################################################

use_col <- c(1:9, which(colnames(seg) %in% info$t_n))
seg_use <- seg[,use_col]
seg_use <- subset( seg_use , Amplitude.Threshold != "Actual Copy Change Given")

###########################################################################################
# Compute the Euclidean distance matrix
dist_matrix <- dist(t(seg_use[,-c(1:9)]), method = "euclidean")

# Perform hierarchical clustering using Ward's method
hclust_result <- hclust(dist_matrix, method = "ward.D")

## 分成两簇
gvhd.cut <- cutree(hclust_result, k = 2)

gs_sample <- names(which(gvhd.cut==1))
cin_sample <- names(which(gvhd.cut==2))

## 
#cin_sample <- cin_sample[!(cin_sample %in% c())]
gs_sample <- gs_sample[!(gs_sample %in% c("JZ585T3_JZ585B" , "S43_S67" , "S52_S70" , "JZGCWES7_JZGCWES0684" , "JZGCWES105_JZGCWES1060" , "JZ661T_JZ661B"))]
## 人工检查分类
## "JZ585T3_JZ585B" , "S43_S67" , "S52_S70", "JZGCWES7_JZGCWES0684"
## 倍体高，同一患者的肿瘤为CIN

## 其判断为GS但是拷贝数改变高
##  "JZGCWES105_JZGCWES1060" , "JZ661T_JZ661B"

#gs_sample <- c(gs_sample , c() )
cin_sample <- c(cin_sample , c("JZ585T3_JZ585B" , "S43_S67" , "S52_S70", "JZGCWES7_JZGCWES0684" , "JZGCWES105_JZGCWES1060" , "JZ661T_JZ661B") )

###########################################################################################

dat_burden$t_n <- paste0( dat_burden$Tumor , "_" , dat_burden$Normal )
dat_burden$CNV_Type <- ""
dat_burden$CNV_Type <- ifelse( dat_burden$t_n %in% gs_sample , "GS" , dat_burden$CNV_Type  )
dat_burden$CNV_Type <- ifelse( dat_burden$t_n %in% cin_sample , "CIN" , dat_burden$CNV_Type  )

###########################################################################################
dat_burden$TCGA_Class <- ""
dat_burden$TCGA_Class <- ifelse( dat_burden$CNV_Type == "GS" , "GS" , dat_burden$TCGA_Class )
dat_burden$TCGA_Class <- ifelse( dat_burden$CNV_Type == "CIN" , "CIN" , dat_burden$TCGA_Class )
dat_burden$TCGA_Class <- ifelse( dat_burden$MS_Type == "MSI" , "MSI" , dat_burden$TCGA_Class )
dat_burden$TCGA_Class <- ifelse( dat_burden$Class == "IM" , "IM" , dat_burden$TCGA_Class )

write.table( dat_burden , out_file , row.names = F , quote = F , sep = "\t" )